Validity of the one-dimensional dissipative Boltzmann equation for point particles up 

to the clustering regime 
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We study stationary states of a one-dimensional gas of point-like particles not subject to gravity 
between two walls at temperatures T_ and T+, with T_ < T+ . Depending on the normalized 
temperature difference A = (T+ — T_)/(T+ + T_) the system may be completely fluidized, or in 
a mixed state in which a cluster coexists with the fluidized gas. We devise and explain in detail 
a method for integrating the one- dimensional dissipative Boltzmann equation in the test-particle 
limit for the stationary case. We then apply this method to test the equation's validity up to 
the clustering regime, by comparing with results from microscopic Newtonian molecular dynamics. 
There is very good agreement, with the one-particle phase space density function presenting highly 
non-Gaussian features, and a discontinuity that corresponds to the test-particle limit. We conclude 
that Boltzmann's equation is valid at least everywhere in the control parameter space where the 
system has no cluster. The behavior of the system in its fluid phase is dominated by characteristic 
lines which resemble trajectories of particles subjected to a force which attracts them to a fixed 
point. If this point is in the physical region a cluster forms, if not then the system remains fluid. 

PACS numbers: 45.70.-n, 05.20.Dd, 02.60.Cb, 02.70.-c 



I. INTRODUCTION 

Granular systems have been the focus of much at- 
tention due both to the theoretical challenges they 
present and to the applications of industrial impor- 
tance that stem from the rich phenomena they exhibit 
(see Refs. JO, @, and references therein). These sys- 
tems are characterized by an energy loss in collisions, 
and this loss is at the base of many interesting phenom- 
ena. Among these, the clustering of particles has drawn 
much attention ||, |, |, @. 

In this paper we study in detail the mechanisms that 
dominate the collective dynamics of a one-dimensional 
system of point-like particles that interact via collisions 
that conserve momentum but dissipate kinetic energy. A 
cluster may or may not form. The system is confined 
in a box of unit length and any particle that reaches a 
wall is expelled from it with its velocity randomly chosen 
so that the velocity distribution of "outgoing" particles 
is a Maxwellian distribution with the "temperature" of 
that wall. There are no external forces. In a previous 
article we saw that there are two relevant control pa- 
rameters: the restitution coefficient which characterizes 
the collisions and the normalized temperature difference 



A between the walls, A = 



T + -T_ 



In the plane of these 



two parameters there is a transition line: on one side 
the system is a granular fluid that reaches a stationary 
regime while in the other side a cluster is formed and, ap- 
parently, no stationary solution can be reached, at least 
in the limit of infinitely many particles. In Ref. we 
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described what happens, while in the present paper we 
disclose the underlying mechanisms. 

Two aspects have to be distinguished. On the one hand 
there is the formal aspect: we show that Boltzmann's 
equation describes closely what we get from molecular 
dynamic simulations in the case of the pure fluid phase 
even quite close to the transition line. On the other, there 
is an intuitive picture — originated in our detailed integra- 
tion of Boltzmann's equation described in this paper — 
that we explain in the following paragraphs. 

To fix notation, if c\ and C2 are the velocities of two 
particles that are about to collide, their velocities after 
the collision are given by 



c[ = qci + (1 - q)c 2 , 



(1 - q)ci + qc 2 . 



Here q = (1 — r)/2, where r is the usual restitution coef- 
ficient. For the elastic case (r — 1) the particles simply 
exchange velocities. Since the grains are point-like, the 
elastic case is then indistinguishable from a system in 
which the particles do not interact. To make this more 
explicit, the point-like character of the grains allows us 
to exchange their identities after the collision, giving the 
collision rules 



qci + (1 - q)c2, 



(1 - q)ci + qc 2 . 
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Thus, when q — the velocities are unaffected, and when 
q is small the velocities are only barely changed. This 
leaves us with the picture of a system of weakly inter- 
acting particles, whose relative velocity diminishes upon 
collisions. 

The one dimensional granular system is being excited 
from the two walls, generally at different temperatures, 
T_ and T+. Particles emerging from the walls act as an 
outcoming "wind" pushing the particles away from them. 
One could picture the effect of this wind as an effective 
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repulsive force which pushes the particles away from the 
walls. If the temperature difference between the two walls 
is large enough, the repulsive force associated to the hot- 
ter wall prevails over the force associated to the colder 
wall all across the system. Therefore in this case the over- 
all effect is a net force always pointing toward the colder 
wall, much like gravity acts in a gas, always pointing to 
the base. If, on the contrary, the temperature difference 
is not large enough, there is a point in the system where 
the two repulsive forces cancel each other, producing an 
equilibrium point — a particle at rest in this point would 
tend to remain at rest — about which a cluster will grow. 
As the cluster absorbs particles the density of the sur- 
rounding gas decreases, and the equilibrium point may 
shift in time. 

The article is organized as follows. In Sec. || we 
state the kinetic equation in the limit which the au- 
thors of Ref. Q call the hydrodynamic limit of the one- 
dimensional Boltzmann equation for point-like grains. 
This limit is analogous to the Boltzmann-Grad limit of 
infinitely many particles and finite mean free path. In 
the present case, the role of the inverse of the mean free 
path is played by the factor qN, where N is the number 
of particles, and q is the inelasticity factor mentioned 
above. This factor represents how much is the velocity of 
a fast test-particle affected by crossing the system, when 
the inelasticity q is very small [Q, ^| . The boundary con- 
ditions and normalization used are also stated. 
Section 



III 



focuses on the stationary state, and de- 
scribes the algorithm devised for the kinetic equation in 
this case. 

Section ^ contains results and conclusions. It com- 
pares the distribution functions in the one-particle phase 
space f(x, c) (where x is the spatial coordinate, and c is 
the velocity) with the same distribution measured from 
molecular dynamic simulations. The numerical solution 
presents a discontinuity that stems from the limit taken 
in Sec. y, which allows the equation to be treated as ap- 
proximately linear. The measured distribution exhibits a 
softened version of the discontinuity that steepens as we 
consider systems with larger number of particles (thus 
approximating better the limit of Sec. |J). An intuitive 
picture is finally put forward. 



II. KINETIC EQUATION, BOUNDARY 
CONDITIONS, AND NORMALIZATION 

In the limit N — > oo, but keeping qN fixed, the 
one-dimensional Boltzmann equation transforms into the 
test- particle equation []|, ^, ||: 



d t f + cd x f = qNd c (Mf), 



(1) 



where 



M (x, c) 



f(x,c')(c-c')\c-c'\dc'. (2) 



Since the equation is nonlinear, we must define explicitly 
the normalization used. In this case it is 



f(x, c) dcdx — 1. 



(3) 



The system is confined in a box of unit length, and the 
particles may have any velocity: (x, c) £ [0, 1] x (— oo, oo). 
Any particle that reaches a wall is instantaneously ex- 
pelled from it — there is no adsorption — with its velocity 
randomly chosen so that the velocity distribution of ex- 
pelled particles is a Gaussian distribution with the tem- 
perature of that wall. This corresponds to choosing a 
wall kernel without memory and without a delay time 
(see Ref. jL(j), for example): 

/(0,c>0)oce- c2 / 2T - /(l,c<0)oce- c2 / 2T +. (4) 

The temperatures at both walls are chosen so that the 
system temperature for the perfectly elastic case is Tq — 
V /T_T+ = I . We will always take T + > T_. 

The missing constants in Eq. (|j) are determined by 
Eq. (H) and by imposing that there is no flow across the 
walls: 



cf(x wa n,c)dc = 0. 



(5) 



III. 



STATIONARY STATE: 
ALGORITHM 



SOLUTION 



In a stationary situation, Eq. ([j]) may be rewritten as 
follows: 



cd x f-qNMd c f = qNfd c M. 



(6) 



The coefficient —qNM multiplying d c f plays the role of 
a force (per unit mass) and it is what we have called 
wind, in the introduction. It is the effective total force 
on a particle at x with velocity c. When / is reasonably 
close to the true solution, M will not depend on the de- 
tailed form of the distribution. Thus, if we have a trial 
distribution /„, we may consider M and dM/dc as given 
functions of x and c, and then we may solve Eq. (^) for 
the distribution f n +\. Seen in this light, when the trial 
function is reasonably close to the solution, Eq. (^) is 
approximately a linear partial differential equation that 
can be analyzed as a hyperbolic equation, using the no- 
tion of characteristic curves jllj . In the present case, the 
characteristic curves satisfy 



dx 

ds 
dc 

— = -qNM(x,c) 
ds 

|J = qNfd c M(x,c), 



(7) 
(8) 
(9) 



where s is a parameter. 
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(a) 




M=0 




FIG. 2: Qualitative picture of the shape of a characteristic 
curve in (x, c, /)-space for the case where the M = curve 
crosses the c = line. 
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FIG. 1: Schematic representation of the form of the char- 
acteristic curves in (x,c)-space when the right wall is hotter 
than the left wall. Subfigure (a) shows the case where the 
curve M(x,c) — does not cross the c = line. Subfigure 
(b) shows how a characteristic curve would wind around the 
point where the M — and c = lines cross. 

In simple words, our integro-diffcrcntial equation is 
treated as if it were a (quasi-linear) partial differential 
equation and, since real characteristics exist, it is pos- 
sible to integrate along these lines as if dealing with an 
ordinary differential equation with a single independent 
variable s. 

Briefly put, given a distribution f n (which implies that 
we have M n and d c M n ), we calculate f n +\ by solving 

c d x f n+ i - qNM n d c f n+1 = qNf n+1 d c M n (10) 

through numerical integration along the characteristics, 
following Ref. jl2|. After the integration we normalize 
fn+i to one, and then use f n +i to calculate f n +2- In this 
way we eventually reach a fixed point. 

As will be described in detail in the following section, 
there are three types of characteristics: (1) those that 
originate at x = with a large positive velocity and end 
at x = 1; (2) those that also begin at x = (implying 
c > 0) but do not reach x = 1: they reach the c = axis, 



turn around, and return to x — 0; (3) those that start at 
x = 1 (implying c < 0) and end at x — 0. The first two 
types of characteristics are associated to the left bound- 
ary condition, while the third type is associated to the 
right boundary condition. The solution, therefore, may 
be discontinuous along the separatrix of these last two 
types of characteristics. Since our numerical algorithm 
integrates along these characteristics, it never crosses the 
discontinuity; every step deals with a smooth function. 

As in Refs. ||, |) we may consider that the projection 
of any characteristic line to the (x, c) plane corresponds 
to the phase-space trajectory of a test particle crossing 
the system. From Eq. (^|) we confirm what has been 
said before, namely that —qNM is the acceleration of 
the particle. For large velocities M ~ c\c\, hence if the 
particle's velocity is large it will be slowed down. This is 
suggested in Fig. |l|(a), where the characteristics far from 
c = approach that axis. 

Due to the form of Eqs. (0) and (||), the intersection of 
the curves c = and M(x,c) = is interesting At 
c = the characteristic curves in the (x, c) plane are ver- 
tical, and at M = they are horizontal. Figure |l](b) 
is a sketch of what would happen to the characteris- 
tic curves around the intersection point (which we will 
call G henceforward): they would wind around it, never 
reaching it. Meanwhile, since d c M > in that vicinity, 
we have that / is increasing along the curve. In other 
words, in [x, c, /)-space the characteristic curve is prac- 
tically vertical, with / increasing sharply along it. This 
case, depicted qualitatively in Fig. ||, corresponds to the 
presence of a cluster at the intersection point G. 

Hence, for the system to be able to reach a stationary 
fluidized state, the curve M = must never cross the c 
axis. In Fig. |l](a) we have drawn this curve in the lower 
part of the plane, which corresponds to the case when the 
left wall is the colder one (T_ < 7+) and the temperature 
difference is large enough that no cluster can form. 

Further details of the algorithm are relayed to the ap- 
pendix. 
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IV. RESULTS AND CONCLUSIONS 




qN 



FIG. 3: Comparison between the threshold for cluster forma- 
tion in molecular dynamics and the loss of convergence in the 
algorithm. The plus signs show the lowest values of A for 
each qN before the algorithm becomes unstable. The circles 
show the lowest values of A for each qN before a cluster is 
detected in a molecular dynamics simulation of N = 1000 
particles. The plus signs are joined by lines to guide the eye. 




FIG. 4: Distribution function f(x, c) for the case qN — 0.35 
and A = 0.6. Only one every four curves is shown to unclutter 
the picture. 

Figure ^ shows the threshold for cluster formation in 
molecular dynamics and the loss of convergence in the 
algorithm. We can say that they coincide, namely, the 
integration of Boltzmann's equation with the present al- 
gorithm converges almost to the transition line beyond 
which a cluster begins to form. 

Figure U shows the characteristic lines in (x, c, /)-space 
when qN = 0.35 and A = 0.6. Although it corresponds 
to a completely fluidized case it is not to far from the 
transition line, as can be seen in Fig. ^. 

Figure shows the projection of the characteristic 
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FIG. 5: Projection of the characteristic curves into (x,c)-space 
for the same case as Fig. [l|. Only one every four curves is 
shown to unclutter the picture. 

curves into the (x,c) plane. The asymmetric aspect of 
the family of curves is due to the effective force which 
pushes particles toward the colder wall almost as if there 
was a space-dependent gravity-like force. The figure only 
shows the characteristic lines for small values of the ve- 
locity since further away the distribution is a bimodal 
Gaussian. One interesting feature of this figure is the 
density of lines near x — for small negative velocity. 
It corresponds to a remarkable peak of the velocity dis- 
tribution for velocities much smaller than the thermal 
velocity. In the present case the thermal velocity for the 
particles approaching the wall is about V2, an order of 
magnitude larger than the width of the peak. 

To check how the distribution obtained from our al- 
gorithm compares with the distribution stemming from 
molecular dynamic simulations we show the distribution 
at the colder wall f(0, c), where the discontinuity is more 
notorious. In Fig. rc| we compare with molecular dynamic 
results for N = 1000 particles, qN = 0.1, and A = 0.6, 
namely quite far from the clustering regime. Figure |^(a) 
shows the distribution for a wide range of velocities (at 
this scale the calculated and the simulated solutions can- 
not be resolved by the naked eye). It can be checked 
that the distribution behaves as two Maxwellians, one 
for c > and a different one for c < (with |c| suffi- 
ciently large), and has a remarkable peak for small neg- 
ative velocities. Even though the system is far from the 
clustering regime, that peak is a reminder that a cluster- 
ing singularity exists. 

Figure ^(b) shows in detail the shape of the discontin- 
uous behavior at the peak. The discontinuity that the 
analytic analysis predicts is softened in the simulations. 
This difference is due to size effects. In fact, as seen 
in Fig. [?] (which corresponds to a system on the verge of 
clusterization) , simulations of systems with an increasing 
number of grains exhibit increasingly steeper behavior at 
the predicted discontinuity, approaching the result that 
Boltzmann's equation implies and our algorithm yields. 
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FIG. 6: Distribution function at the cold wall for the case 
qN — 0.1 and A = 0.6. The solid curve results from ap- 
plying the algorithm described in the text, while the dashed 
curve was measured from a molecular dynamics simulation 
with N — 1000 particles. Subfigures (a) and (b) show the 
same data, but on a different scale. Subfigure (a) empha- 
sizes how the distribution is essentially Gaussian, save for the 
sharp peak for slow particles. Subfigure (b) shows a detail of 
this peak, showing how the discontinuity is smoothed in the 
simulation. 



The final picture that emerges from all that has been 
said is, first, that Boltzmann's equation for the quasi- 
elastic system is valid essentially in all points of the (q, A) 
plane where the system has no cluster. The behavior of 
the system in its fluid phase is dominated by character- 
istic lines, trajectories of a test particle subjected to a 
force which attracts it to a point G where the M = 
line crosses the c-axis. In the fluid phase such point G 
is beyond the physical box, particles hit the colder wall, 
forget their past and reenter the system. 

We have not solved the time-dependent case when G 
is inside the box. In such case many trajectories in phase 
space will wind around G (as shown in Figs. l](b) and|^), 
and the density at that point will tend to diverge, thus 
forming a cluster. 
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FIG. 7: Detail of the discontinuity in the distribution function 
at the left wall for qN = 0.35 and A = 0.6. The thick solid line 
results of applying the algorithm described in the text. The 
other curves stem from molecular dynamics for 1000 (dashed) , 
2000 (dash-dotted), 5000 (dotted), and 10000 (thin solid line) 
particles. 
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APPENDIX: ALGORITHM IN DETAIL 

In order to commence the integration, we need some 
initial ansatz fo(x,c). This may be, for example, the 
solution for the elastic case: 



f(x, 



A 



{e(c) ( 



2 /2T_ 



B[l - 9{c)]c 



2 1 IT. 



}• 



where 9(c) is the Heaviside step function, B is a constant 
chosen so as to have zero flux at the walls, and A is a 
normalization constant. In practice, if it is available, it is 
convenient to choose as fo a previous solution for a case 
similar to the one being studied. This not only speeds 
up the convergence, but also may help avoid the spurious 
appearance of G points at intermediate iterations. 

The next step is to see whether the M = curve crosses 
the c = line. If it does, the algorithm breaks down, and 
we must start from a better ansatz. If it lies wholly on 
c < 0, the characteristic curves will behave qualitatively 
as in Fig. |l](a). In this case we start by integrating the 
characteristics that begin at (x = 0, c > 0). If the curve 
lies wholly on c > the situation is inverted, and we must 
start from (x — 1, c < 0). In what follows we will assume 
that the situation is as in Fig. 0(a). This is always the 
case when we are sufficiently close to the solution. 

The next step is finding the dividing characteristic. 
This is done directly by the bisection method, choosing 
different values of cd so that the initial condition for the 
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characteristic is 

x(s = 0) = 0, c(0) = c D , /i(0) = A-e- c °' 2T - . 

The value of A- is chosen so that the net flux at the wall 
is zero: 

/0 />oc 
cf Q (0,c)dc + A- / ce- c2/2T -dc = 0. 
-oo J0 

Integrating along this characteristic in very small steps 
we find cd such that the corresponding curve has a turn- 
ing point at (a;, c) = (1 — e, 0), with e less than a reason- 
ably small number. 

With this value of cjj we divide the interval < c < cjj , 
choosing many values c, (with < ci < C2 < . . . < cd) 
as starting points for characteristic curves. These will all 
be curves which will turn around and return to a; = 0. 
The crossing points will define a natural discretization of 
the x axis which will be used to tabulate the values of 
fi. In this way we always sample the turning point. 

Now we start the integration: with increasing i, we in- 
tegrate the characteristic that starts with Cj, tabulating 
the values of c(s) and f(s) for values of x corresponding 
to the turning points of the previously integrated char- 
acteristics. This step is of paramount computational im- 
portance: at the next iteration we will need to calculate 
the integrals M(x,c) and d c M(x,c). Since these are in- 
tegrals in c keeping x fixed, it is important to have values 



of f(x,c) tabulated at the same values of x. This pre- 
caution allows us to evaluate M and d c M in a fast and 
straightforward manner. 

Having integrated along all the characteristics with 
Ci < cd, we now integrate characteristics that start from 
(x = 0, c > Cd)- These are curves that start at the left 
wall and reach the right wall, that is, they do not have a 
turning point. The completion of this step implies that 
we have /i in the whole half-space c > (and also in part 
of c < 0)'. 

At this point we can start integrating the characteristic 
curves that start at the right wall, choosing values of 
c < 0: 

x(s = 0) = l, c(0) = c , .M0) = A+e- c o/ 2T +. 

The value of A + is chosen so that the net flux at the wall 
is zero: 

/0 />oc 
ce- c2/2T +dc+ / c/i(l,c)dc = 0. 
-oo JO 

After completing the integration, we normalize /i, and 
take the magnitude of this adjustment as a measure of 
how far we are from reaching a fixed point. 

This procedure is repeated until a reasonable conver- 
gence is reached. Typically ten to fifteen iterations are 
necessary to achieve a good solution. 
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